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ABSTRACT 

We have used the hydrodynamical AMR code ENZO to investigate the dynamical 
evolution of the gas at the centre of dark matter haloes with virial velocities of ~ 
20 — 30kms _1 and virial temperatures of ~ 13000 — 300007^ at z ~ 15 in a cosmological 
context. The virial temperature of the dark matter haloes is above the threshold where 
atomic cooling by hydrogen allows the gas to cool and collapse. We neglect cooling 
by molecular hydrogen and metals, as may be plausible if H2 cooling is suppressed 
by a meta-galactic Lyman- Werner background or an internal source of Lyman- Werner 
photons, and metal enrichment has not progressed very far. The gas in the haloes 
becomes gravitationally unstable and develops turbulent velocities comparable to the 
virial velocities of the dark matter haloes. Within a few dynamical times it settles 
into a nearly isothermal density profile over many decades in radius losing most of its 
angular momentum in the process. About 0.1 - 1% of the baryons, at the centre of 
the dark matter haloes, collapse into a self-gravitating, fat, ellipsoidal, centrifugally 
supported exponential disc with scale-length of ~ 0.075—0.27 pc and rotation velocities 
of 25 — 60 kms -1 . We are able to follow the settling of the gas into centrifugal support 
and the dynamical evolution of the compact disc in each dark matter halo for a few 
dynamical times. The dynamical evolution of the gas at the centre of the haloes is 
complex. In one of the haloes the gas at the centre fragments into a triple system 
leading to strong tidal perturbations and eventually to the in-fall of a secondary smaller 
clump into the most massive primary clump. The formation of centrifugally supported 
self-gravitating massive discs is likely to be an important intermediary stage en route 
to the formation of a massive black hole seed. 

Key words: Cosmology: theory - large-scale structure - black holes physics - meth- 
ods: numerical 



1 INTRODUCTION 

Supermassive black holes (SMBH)s were invoked as the 
central engine powering quasi-stella r objects (QSO)s soon 
after the first QSOs were discovered ([ZerDovich &; Novikovl 
119641 : ISalpete3ll964h . 

Early predictions that supermassive 
black holes are ubiquitous and maybe found in many if not 
most galaxies as the remnants of dead QSOs ([Lvnden- Belli 
Il969h have stood the test of time. SMBHs are now believed 
to reside in most if not all galactic bulges and their masses 
correlate tightly with the stellar velocity dispersion and 
the mass of the bulges of their host galaxies (Gebhard 
et al. 2000; iFerrarese fe Merrittl l2000h . The luminosity 
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functions of active galactic nuclei across the electromag- 
netic spectrum and its evolution with redshift have been 
used to constrai n the growth history of supermassive 
black holes (e.g.. IMerloni et al.1 120041 : IShankar et all 120071 : 
iMerloni &; Heinzll2008h . Little is known however, about how 
the first (super-) massive black holes came into being. The 
discovery of very luminous, bright QSOs at z > 6 appears 
to suggest that at least some of the most massive black 
holes were already in place when the Universe was less 
than 1 billion years old |Fan 2001; Fan et alJl20o3 : iHaimanl 
l2Q06h . 

A wide range of generic pathways to a supermassive 
black hole are possible ranging from direct collapse of 
gas into rather massive seed black holes, to Eddington- 
limited accretion onto stellar mass black holes and the 
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Sim 


Boxsize 




Z co \\ 


DM mass 


A R 






[Comoving h 1 Mpc] 






[M ] 


[PC] 




A 


5.0 


200.0 


15.0 


9.58 x 10 3 


9.59 x 10 -3 




B 


2.0 


250.0 


14.0 


4.42 x 10 2 


1.61 x 10 -2 




C 


2.0 


250.0 


15.0 


4.42 x 10 2 


1.61 x 10~ 2 




Mtot 


R-200 


V c 


T v ir 


Pmax 


A 


Tcore 


[M ] 


[kpc] 


[km s- 1 ] 


[K] 


[cm- 3 ] 




[K] 


2.64 x 10 8 


1.28 


29.9 


3.23 x 10 4 


1.24 x 10 9 


0.031 


6.11 x 10 3 


5.37 x 10 7 


0.64 


19.0 


1.30 x 10 4 


5.86 x 10 8 


0.026 


6.35 x 10 3 


5.15 x 10 7 


0.59 


19.3 


1.35 x 10 4 


7.12 x 10 8 


0.050 


6.13 x 10 3 



Table 1. Basic properties of the three simulations (A, B and C): boxsize (comoving h~ 1 kpc), starting redshift, collapse redshift, DM 
particle mass (Mq), spatial resolution (pc), total mass of the halo (Mq), virial radius (h~ 1 kpc), circular velocity (km s — 1 ), virial 
temperature (K), maximum baryon density in the halo (cm -3 ), angular momentum parameter A, temperature at the core of the halo 
(K). All units are physical units, unless explicitly stated otherwise. 



dyna mical evolution of dense star clusters (|Reesl Il978l 
Il984h . The tight relation between galactic bulges and 
their central mass suggests that they have formed and 
grown in a connected way. In the well-established ACDM 
paradigm of structure formation, dark matter haloes 
grow by hierarchical merging. The last decade has seen 
extensive and detailed modeling of the build-up of galax- 
ies and supermassive black holes in this context using 



(lEfstathiou & Reeslll988l: ICarlberd Il99d 


; lHaehnelt & Reed 


19931: iKauffmann k Haehnelt 2000 


; iMonaco et al.1 


20001: Cavaliere k Vittorinil 120021: IVolonteri et al. 120031: 


Volonteri & Pernal 120051: IVolonteri & Reesll2005l: Volonteri 



2008; lHopkins et al.ll2008f ). Due to the lack of observational 
constraints it is, however, still an open question at what 
mass scale this hierarchical build-up starts. 

The modeling of the growth of supermassive black holes 
by hierarchical merging of dark matter haloes/galaxies 
from the remnant black holes of the first generation of 
stars has identified a number of problems. The shallow 
potentials hosting the first generation of sta rs are rather 
sensitive to the effect of stellar feedback (iDekel Silk! 
Il986l : IWise fe Abel l2007al : lO'Shea fe Normanl l2008h and 
may thus not provide sufficient amounts of fuel for a rapid 
growth of stellar mass black holes. The predicted frequent 
merging may lead to the expulsion due to gravitational 
recoil and/or gravitational slingshot unless stel l ar see d 
black holes are sufficiently rare (e.g. IVolonteril l2Q07bh . 
The available time for the formation of the most massive 
black holes requires almost continuous accretion at the 
Eddington rate or super- Eddington accretion rates. There 
is also circumstantial evidence that a certain class of X-ray 
sources, Ultra- luminous X-ray sources or ULXs for short, 
are powered by i ntermediate mass black holes of ma sses 
up t o 10 4 (pvlakishima et all l200d: IFabbiano et "all 

120031 : IColbert et all 12004 IFabbiano et al.1 l20Q6h . There 
is thus plenty of motivation to consider seriously the 
possibility that the growth of supermassive black holes 
has started from massive seed black holes with masses 
subs tantially larger than t hose forming as stellar remnants 
(see iBegelman et al.l (|2006h for a recent discussion) . 

The centres of high-redshift dark matter haloes 



with virial temperatures > 10 4 K, not yet significantly 
enriched with metals, have been identified as a promis- 
ing environment to form such massive seed black holes. 
If H2 cooling is suppressed in these haloes either due 
to an external UV background or more likely espe- 
cially in the later stages of t he collapse due to i nternal 
sources of UV radiation (cf. IWise fe Abel l20Q7bh early 
fragmentation should not occur favouring the forma- 



(lOh & Haiman 120021: iBromm & ] 


Loebll2003l: IVolonteri et al. 


20031: iKoushiaoDas et al.1 12004 


; Begelman et al.l 


2006; 


Lodato & Nataraianl 120061: 


Rees & Volonteril 


2007; 


Volonteri et al.l l2008h . This is rather different from 



the situation in the lower mass haloes studied extensively 
as possible sites for the formation of the "first" stars 
where H2 cooling is th e dom i nant cooling mechanism 



g. IBromm et al.l Il999l 120021 : 



IBromm fc Larso 
Il998l . I2000L 120021 : lO'Shea &; Normal 



2004 



2007, 



Abel et al 
20081 ). 

We use the publicly available code EnzcQ to perform 
adaptive mesh simulations of the dynamical evolution of 
the gas in three dark matter haloes with virial velocities 
between ~ 19 kms -1 and ~ 30 kms -1 and virial temper- 
atures between ~ 13000 K and ~ 30000 K at z ~ 15 in 
a cosmologic al con t ext. T he work is similar in spirit to 
that of IWise et al.1 <l2Q07h , hereafter W07, who simulated 
two haloes somewhat less massive than our haloes. W07 
find that the gas at the centre of their haloes does not 
reach rotational support and does not fragment. Our 
simulations have a similar set-up but unlike W07 we do 
not push for resolution in the central region where a very 
small fraction of the gas can collapse to very high density 
but rather follow the dynamical evolution of a significant 
fraction of the gas settling into a rotationally supported 
disc for several dynamical times. The paper is structured 
as follows. In 32] we describe the details of the numerical 
simulations. In ^3]we summarise some basic formulae that 
we use to characterise the properties of the haloes. In 2] 
we describe the results of our numerical simulations and in 
35] we summarise our conclusions. Throughout this paper 
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we assume a standard ACDM co smology with the following 
parameters (jSpergel et all I2QQ3I . based on the WMAP 1st 
year data), Q A , = 0.74, Q m , = 0.26, Q b = 0.0463, a 8 = 
0.9 and h = 0.72. We further assume a spectral index of 
primordial density fluctuations of n = 1. 



2 THE SETUP OF THE NUMERICAL 
SIMULATIONS 

2.1 The adaptive- mesh refinement code Enzo 

We have used the publicly available adaptive mesh refine- 
ment (AMR) code Enzo. Enzo was originally developed by 
Greg Bryan and Mike Norman at the University of Illinois 
(Bryan & Norman 1995b, Bryan & Norman 1997, Norman 
& Bryan 1999, O'Shea et al. 2004). The gravity solver in 
Enzo uses an N-Body particle mesh technique (Efstathiou 
et al. 1985, Hockney & Eastwood 1988). Additional finer 
meshes can then be added in regions of high density to cal- 
culate the dynamics of the dark matter (DM) particles more 
accurately. 

The hydrodynamics solver employs the Piecewise Parabolic 
method combined with a non-linear Riemann solver for 
shock captu ring. The Eulerian AM R scheme was first de- 
veloped by Berger &; Oligerl (| 19841 ) and later refined by 
iBerger fe Colellal (1989[) to solve the hydrodynamical equa- 
tions for an ideal gas. Bryan & Norman (1997) adopted such 
a scheme for cosmological simulations. In addition to this 
there are also modules available which compute the radiative 
cooling of the gas together with a multispecies chemical reac- 
tion network. There are two versions of the chemistry solver 
available, one with 6 species (H, H + , He, He + , He ++ , e~) and 
one with 9 species (same as before plus H2,H^,H _ ). As 
stated previously the simulations conducted here do not 
include H2 cooling and chemistry. Our simulations make 
extensive use of Enzo's capability to employ nested grids 
which we describe in the next section. 



2.2 Nested grids &; initial conditions 

Initial conditions were generated with the initial conditions 
generator supplied with the Enzo code. The nested grids 
are introduced at the initial conditions stage. We have first 
run exploratory DM only simulations with coarse resolu- 
tion, setting the maximum refinement level to 4. These DM 
only simulations have a root grid size of 128 3 and no nested 
grids. In these exploratory simulations we have identified 
the most massive halo at a redshift of 10 and then rerun the 
simulations, including the hydrodynamics module. We also 
introduce nested grids at this point. The nested grids are 
placed around the region of interest, as identified from the 
coarse DM simulation. We have used four levels of nested 
grids in our simulations with a maximum effective resolu- 
tion of 1024 3 . The introduction of nested grids is accom- 
panied by a corresponding increase in the DM resolution 
by increasing the number of particles in the region of in- 
terest. We only use the highest resolution nested grid to 
refine further thus economising our computational output 
.This corresponds to a region of 625 h~ x (comoving) kpc for 
simulation A and 250 hT 1 (comoving) kpc for simulations 
B & C. The total number of particles in our simulation is 



4935680, with 128 3 of these in our highest resolution region. 
The nested grids are distributed as follows in root grid sizes, 
L0 = 128 3 , LI = 32 3 , L2 = 24 3 , L3 = 16 3 . Table [Q gives the 
details of the simulations discussed here. 



2.3 Following the collapse of the gas and dark 
matter with Enzo 

As mentioned previously Enzo uses adaptive grids to pro- 
vide increased resolution where it is required. For the simu- 
lations discussed in this paper we have used four refinement 
criteria implemented in ENZO: DM over-density, baryon over- 
density, Jeans length and cooling time. The first two crite- 
ria introduce additional meshes when the over-density of a 
grid cell with respect to the mean density exceeds 3.0 for 
baryon s and/or DM. The th ird criterion is the Truelove cri- 
terion (|Truelove et al.lll997h which in its original form states 
that at least 4 grid cells should be used to resolve the Jeans 
length to ensure that no artificial fragmentation takes place. 
The cooling time is often shorter than the dynamical time 
and the cooling time refinement criterion helps to ensure 
that the Jeans mass is prope rly resolved. As other authors 
(e.g. lO 'Shea &; NormanlEooil . W07) we are conservative here 
and set this criterion to 16 grid cells in our simulations. The 
fourth criterion ensures that the cooling time in a given grid 
cell is always longer than the sound crossing time of that 
cell. We also set the MinimumM ass For Refinement Exponent 
parameter to —0.2 making the simulation super-Lagrangian 
and therefore reducin g the criteria for refinem ent as higher 
densities are reached (jO'Shea fc Normanll2008n . We further- 
more set the MinimumPressureSupportP ammeter equal to 
ten as we have rest ricted the maximum ref inement level in 
our simulations fe.g. lKuhlen &; Madaull2005h . With the Min- 
imumPressureSupport option the code assumes in each com- 
putational cell the minimum temperature necessary to make 
the cell Jeans stable at the highest refinement level. 

We first advance our simulation forward in time using 
a maximum refinement level of 4. If we would have run the 
simulation with a high maximum refinement level from the 
outset the simulation would have stalled once the virial tem- 
perature of the halo has reached 10000K in the attempt 
to follow the dynamical evolution of a small dense region 
before a sufficiently massive DM halo had formed. We are 
effectively delaying the collapse by limiting the resolution. 
We use a friends of friends (FoF) algorithm to keep track of 
how our halo is progressing through time. Once a sufficiently 
massive halo (corresponding to a DM particle number of 
about 30000 particles in simulation A and 120000 particles 
in simulations B and C) has formed we halt the simulation. 
We then restart the simulation with the maximum level of 
refinement set to 18 (16 in simulation B and C). We will 
refer to this point as the collapse redshift. 

The simulation then continues at a much higher reso- 
lution albeit at a significantly slower speed. We allow the 
simulation to evolve until the code no longer tracks the hy- 
drodynamic evolution of the gas accurately due to a lack 
of resolution at the highest densities caused by restricting 
the highest allowed refinement level. The Enzo code issues 
warnings alerting the user when the code begins to produce 
spurious results due to a lack of resolution. At this point we 
terminate the simulation. 

We have also followed the recommendation in 
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Figure 1. A sequence of visualisations of the density distribution in simulation A. The gas is collapsing in a halo of DM mass of 
2.64 x 10 8 Mq, virial velocity of ~ 30kms — 1 and virial temperature of ~ 32000 K. At the end of the simulation the gas at the centre of 
the halo has settled into a rotationally supported disc. Each plot shows a thin slice of the density distribution centered on the grid cell 
with the highest density in our halo of interest. Between panels either the density range or the time of the simulation output change. 
The colour range represents approximately an order of magnitude range in density in each plot. The scales are proper distances. 
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Time = 6.908 Myrs 



Figure 2. Same as figure [T]for simulation B. The gas is collapsing in a halo of DM mass of 5.37 x 1O 7 M0, virial velocity of ~ 19kms _1 
and virial temperature of ~ 13000 K. Note how the gas fragments into three clumps at T = 6.442 Myrs which tidally distort each other 
and undergo a violent dynamical interaction. 
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Time - 50.103 Myrs 



2 pc 



Time - 50.274 Myrs 



Time - 50.550 Myrs 



Time - 50.481 Myrs 




Figure 3. Same as figure [T] for simulation C but only for the later stages of the evolution. The gas is collapsing in a halo of DM mass 
of 5.15 x 1O 7 M0, virial velocity of ~ 19.3 kms - 1 and virial temperature of ~ 13500 K. Panels 9 to 15 illustrate the dynamical evolution 
of the rotationally supported gas at the centre of the halo for a few dynamical times. 
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Figure 4. Left-hand Panels: The temperature profile for the halo in simulations A and B as a function of radius. Middle Panels: The 
time evolution of the radial velocity (black), thermal velocity (blue) and the turbulent velocity (red) as a function of total enclosed 
gas mass. Right-hand Panels: Ratio of the gravitational acceleration to the acceleration due to thermal pressure as a function of total 
enclosed mass. Note that initially the central regions are supported by thermal pressure. 



iLukic et al.l (|2007l ) for the starting redshift of our simula- 
tions. As a result the initial redshift is quite a bit higher 
than that of other simulations described in the literature. 
Experimenting with the value of the initial redshift suggests, 
however, that this increased initial redshift has little effect 
on our results. 



3 CHARACTERISTIC FORMULAE 

As mentioned before we use a FoF algorithm to identify 
DM haloes in our simulation outputs. We adopt a link- 
ing l e ngth of 0.2 for al l of our analysis (see Ijenkins et al.l 
l200ll ; | Lukic et al. I l2Q07h . We compute physical properties 
of the virialised h aloes using standard formulae (see e.g. 
iMo Whitell2002h . For convenience we summarise some of 
the formulae below. We calculate characteristic properties 
for a sphere for which the mean enclosed density is 200 times 
the mean cosmic value ~p. The characteristic "virial" radius 
is then related to the mass as, 



R200 



GM 



100Qm(z) H 2 (z)_ 
while the circular velocity, V c , can be written as 

, GM(r) 



(1) 



(2) 



where H(z) is Hubble's constant at redshift z and Q m (^) is 
the density parameter of non-relativist ic matter, M is the 
mass of our halo as determined by FoF group finder. Rewrit- 
ing the equation for the virial velocity in terms of the mass, 
M, only gives 



V c = M 1/3 ( y/l00n m (z)GH(z) 



1/3 



(3) 



The Hubble constant and density parameter are related 
to their present-day values by 



H{z) = HoE(z) 
and 



(4) 
(5) 



EP(z) ' 
where E(z) is given by 

E(z) = [QA,o + (l-n m ,o-nA,o)(l + ^) 2 + n m ,o(l + ^) 3 ] 3 -(6) 
We can now define the virial temperature of the halo as 

Vc 



2k 



3.24 x 10" 



30 kms" 



K, 



(7) 



where \i — 0.6m p , with m p being the proton mass, is the 
mean molecular weight. 

The level of rotational support is often characterized 
by the dimensionless a ngular momentum p arameter A which 
can be written as (e.g. iBullock et alJlioQlh 



A 



y/2MV c r ' 



(8) 



where | L | is the angular momentum inside a sphere of radius 
r containing mass M and V c is the virial velocity of the halo. 
The mean value of A is ^ 0.04 for typical haloes in cosmo- 
logica l simulations ([Barnes &; Efstathioul Il987l ; iBett et all 
120071 ). We have included the above properties for the three 



© 0000 RAS, MNRAS 000, 000-000 



8 J. A. Regan et al 





Figure 5. Left-hand Panels: The evolution of the enclosed (gas) mass within spherical shells around the centre of the halo as a function 
of radius. The enclosed DM mass profile at the end of the simulation is shown by the filled diamonds. Right-hand Panels: The evolution 
of the (gas) density profile during the collapse. The gas settles into an close to isothermal (r -2 ) density profile (thick dot-dashed line) 
over many decades in radius. The DM density at the end of the simulation is shown by the filled diamonds. 



virialised haloes in table [T] for reference. 

Also of interest is the dynamical time which we express 
in terms of the enclosed mass M(r) at radius r, 



4 RESULTS OF THE NUMERICAL 
SIMULATION 

4.1 The global dynamics of the collapse of the gas 
and dark matter 

We first discuss the global dynamical evolution of the gas 
and dark matter in the three haloes we have simulated. Note 
that the dynamical evolution in simulation C is similar to 



that in simulation A and we will thus show only a sub- 
set of the plots for simulation C in the following. We have 
simulated boxes with 5/i _1 Mpc and 2/i _1 Mpc on a side, 
respectively. Figure 1 and 2 show the location of the haloes 
in simulation A and B within the typical characteristic web 
of filaments and sheets which is already well developed at 
z ~ 15 (simulation C looks very similar in this respect). 
The haloes are located at the intersection of several fila- 
ments. In the top row of figure 1 and 2 we zoom in several 
steps onto the haloes which are approximately spherical and 
have a virial radius of about 1.5 and 0.5 kpc. In Figure 3 
we focus on the later stages of the evolution of the gas at 
the centre of simulation C. From one panel to the next we 
either change the output time or the density /length scale. 
We have selected time as the time at which we re-start 
the simulation with increased resolution, with a maximum 
refinement level of 18 for simulation A and a maximum re- 
finement level of 16 for simulations B & C. Recall that up 
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Figure 6. : Left-hand Panel: The evolution of the total specific angular momentum (I) of the gas as a function of enclosed mass for 
eight of the simulation outputs in figure^ Right-hand Panel: The same for simulation B for eight simulation outputs shown in figure 

m 



to now we had run the simulation limiting the refinement 
to 4 levels in order to allow the halo to build up. We will 
discuss possible limitations of this approach later. Once the 
resolution is increased the inner part of the haloes starts to 
collapse. 

As discussed by W07 the collapse is highly turbulent 
and dynamically complex. In each halo we are able to fol- 
low the collapse of the gravitationally unstable gas at the 
centre of the halo into a centrifugally supported disc (or 
disc-like object) and follow the dynamical evolution of the 
disc for several dynamical times. In simulation A this disc 
has a mass of 2 x 10 4 M©. In the second of our haloes the gas 
at the centre of the halo fragments into three gravitationally 
bound clumps with masses of a few times 10 4 M© and a mass 
ratio of approximately 3:1:1. The clumps tidally distort each 
other and subsequently undergo a violent dynamical inter- 
action. One of the smaller clumps eventually merges with 
the most massive clump and forms a disc of ~ 1 x 1O 5 M0. 
The second small clump is tidally stripped of most of its 
mass and has a mass of a few times 10 3 M© at the end of the 
simulation. The (temporary) formation of multiple systems 
appears, however, not to be generic. Neither simulation A 
with its five times more massive halo nor simulation C with 
its halo of similar mass show the formation of a similar mul- 
tiple system. The disc forming in simulation C has a mass 
of - 1 x 1O 5 M . 



4.2 The onset of gravitational instability, mass 

inflow and the evolution of the density profile 

The gas at the centre of our haloes attains an almost isother- 
mal temperature profile with a characteristic temperature 
of 6000-7000K the temperature below which atomic cooling 
due to hydrogen cooling becomes inefficient. The tempera- 
ture rises gently to 9000-10000 K at larger radii. The left- 
hand panels in figure |4] show the temperature profile of the 



gas in the haloes and its evolution for simulation A (top 
panel) & B (bottom panel), respectively. 
In order to understand the dynamical evolution of the gas it 
is illustrative to consider the thermal and turbulent veloci- 
ties of the gas. In the middle panel of figure |4] we show the 
evolution of the radial velocity (black) , the thermal velocity 
(blue) and the turbulent velocity (red) of the gas. As the 
collapse gets under way the gas develops turbulent veloci- 
ties comparable to the virial velocities of the halo starting 
in the outer parts of the halo and progressively moving in- 
wards. The turbulent velocities are calculated by computing 
the root mean square velocity of the gas after subtracting 
the centre of mass velocity of the halo and the velocity due 
to the radial inflow of the gas. We also subtract an estimate 
of the rotation velocity (based on the specific angular mo- 
mentum and the inertia tensor as described in section 4.3) 
in quadrature. The turbulent velocities in the later stages of 
the collapse become significantly larger than the thermal ve- 
locities of the gas which we calculate as Vth = \/SkT/ firriH , 
where fi is the mean molecular weight, chosen to be 1.22 and 
rriH is the mass of the hydrogen atom. As in the simulations 
of W07 the turbulent velocities are supersonic. There is a 
net inflow of gas with velocities about a factor three smaller 
than the turbulent velocity component. 

In the right hand panel of figure |4] we compare the ra- 
tio of the inward acceleration due to gravity to the out- 
ward acceleration due to the thermal pressure (gradient) for 
both simulation A (upper panel) and simulation B (bottom 
panel). Once the DM halo has built up, the mass of the 
halo is above the Jeans mass and the inward gravitational 
acceleration comfortably exceeds the outward acceleration 
due to the thermal pressure. The peak which builds up at 
an enclosed mass of a few times 10 4 Mq for simulation A 
and ~ 10 5 Mq for simulation B is due to the formation of a 
centrifugally supported disc at the centre. The double peak 
at certain output times for simulation B is due to the pres- 
ence of multiple gas clumps within the system. 
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Figure 7. Left-hand Panels: The circular velocity V^(blue), and our two estimates of the rotation velocity / /a\ (black) and I /r 
(red) are plotted as a function of radius. As the collapse proceeds the inner part of the halo settles into centrifugal support. Right-hand 
Panels: The ratio of V ro t = I /cl\ to V c , for the same six output times plotted against radius. The central region attains rotational 
support, after w 12.6 Myrs and « 6.7 Myrs, in simulations A and B, respectively. 



In figure [5] we show the evolution of the enclosed gas 
and DM mass as a function of radius and the density pro- 
files of the gas and DM distribution. The DM density (in 
units of ra # cm -3 ) and the DM fraction are shown for the 
final output times only as filled diamonds. The density pro- 
files are calculated by averaging over spherical shells centred 
on the densest point in the halo. Note that in simulation B 
the densest point is always found in the most massive of the 
three clumps. The haloes collapse within 12 and 7 Myrs, re- 
spectively, after we have restarted the simulation with the 
increased refinement level. This corresponds to about a few 
free-fall times for the gas in the inner few hundred parsecs 
which initially has pretty much constant density. As dis- 
cussed already, we initially held the refinement level at a 
maximum of four until the halo had a mass corresponding 
to - 30000 ( - 120000) DM particles. The gas completely 
decouples from the dark matter during the collapse and be- 
comes self- gravitating. The density profile of the gas in the 



inner part steepens dramatically and settles into a close to 
r~ 2 density profile over many decades in rad ius as expected 
for an isothermal collapse (e.g. lLarsonlll969l W07). The in- 
ner few times 10 4 M (1O 5 M ) collapse further and settle 
eventually into rotational support. In the bottom right hand 
panel of figure [5] (simulation B) the secondary peak in the 
density profile is due to the secondary clump which even- 
tually merges with the primary clump. The peak is absent 
in the final two output times as one of the smaller clumps 
has merged with the primary clump and the other smaller 
clump has had most of its mass stripped away. Note that 
the continuous accretion onto the virialised haloes from the 
filaments occurs on much longer time scales and thus has no 
visible effect. 
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Figure 8. The ratio of the square root of the smallest eigenvalue, 0,3, of the inertia tensor to the distance to the centre of the halo. 
Note the dip at a few times 10 4 Mq where the gas settles into a rotationally supported fat disc in simulation A. For simulation B the 
dip at the end of the simulation occurs at a much higher mass of ~ 1 x 10 5 Mq. 



4.3 Angular momentum loss and settling into 
rotational support 

In Figure [6] we show the evolution of the specific angular 
momentum, |/|, as a function of enclosed mass for both sim- 
ulation A and simulation B. The angular momentum vector 
is calculated in the usual way via the cross product L = fxp 
where p is the momentum vector. We then obtain the spe- 
cific angular momentum vector as a function of enclosed 
mass centred on the centre of mass by dividing by the en- 
closed mass. In both haloes there is very significant angular 
momentum loss. In simulation A the angular momentum 
drops by a factor of 20 or more between the initial output 
and the final output in the inner few times 1O 4 M . The 
evolution in simulation B is more complicated due to the 
complex dynamical interaction of the triple system. Initially 
the angular momentum of the gas drops by a factor 20-100 
in the primary clump and then increases again by a factor 
three to five when one of the smaller clumps merges with 
the primary clump. 

We now want to discuss in more detail to what extent 
the gas settles into rotational support. During the turbulent 
collapse of the gas, it is not obvious how best to define ro- 
tation velocities. We follow W0 7 who use the ratio 
where / is the specific angular momentum vector and r is 
the position vector from the centre of mass as a simple but 
rough measure of the rotation velocity of the gas. We also 
calculate a second estimate of the rotation velocity based on 
the angular momentum / and the inertia tensor, /. 

The nine components of the inertia tensor are given by 

Ixx = ^2rnj(yj +z*), (10) 
3 

Ixy = - ^mjXjVj, (11) 
3 

and the corresponding cyclic permutations, where the 
sum is over computational cells. The off-diagonal compo- 



nents are called the products of inertia while the diagonal 
components are referred to as the moments of inertia. The 
matrix is symmetric which guarantees that the eigenvalues 
are real. 

The angular momentum and the inertia tensor are re- 
lated as, 

1 = Iu, (12) 

where cj is the angular velocity. Using the square root 
of the largest eigenvalue of the inertia tensor, ai, we then 
estimate the rotation velocity as 

Kot^^H. (13) 

In figure [7] we compare our two estimates of the rotation ve- 
locity \l\l\r\ (red) and \l\/a± (black) to the circular velocity 
(blue), V c (r) = sj GM(r)/r, for simulation A (top panel) 
and simulation B (bottom panel). Early on in the collapse 
the gravitational potential is dominated by the dark mat- 
ter halo and the gas is only slowly rotating with a ratio 
of rotation to circular velocity of about 1 : 20. As the gas 
in the inner part of the haloes collapses and becomes self- 
gravitating both the circular velocities and the rotation ve- 
locities rise. As shown in the right hand panel, the inner few 
times 10 4 Mq reach (approximate) rotational support with 
\l\/ai slightly larger than V c after about 12 Myrs in sim- 
ulation A. In simulation B the increase in specific angular 
momentum at around 6.5 Myrs due to the merger of one of 
the small clumps with the primary clump leads to a sharp 
increase in the rotation velocity. As we will see in more detail 
later the gas nevertheless settles into a regular rotationally 
supported disc despite this rather violent dynamical evo- 
lution. Our estimate of the rotation velocity based on the 
inertia tensor exceeds the circular velocity by a factor of up 
to two. We will come back to this point in section 4.5. 

Note that our two estimates of the rotation velocity 
trace each other, but that \l\/\r\ is systematically lower than 
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Figure 9. A close up of the final panels of Figures 1 and 2 showing the 2D density projection of the rotationally supported gas at the 
centre of haloes at the end of simulation A and B. 



by a factor 1.5-3. As the gas settles into rotational sup- 
port a disc-like flattened structure should gradually form. 
For a turbulent collapse like the one considered here we 
found this to be most easily studied by looking at the evo- 
lution of the smallest eigenvalue of the inertia tensor which 
for a disc should be a good proxy for the thickness of the 
disc. In figure [8] we show the ratio of the square root of the 
smallest eigenvalue of the inertia tensor, <23, to the radius 
in which the inertia tensor is calculated as a function of the 
enclosed mass. The ratio reaches a minimum value at a few 
times 10 4 M© and a few times 10 5 M© for simulation A and 
B, respectively, as expected for the formation of a rotation- 
ally supported disc. The disc appears to fatten towards the 
centre and is surrounded by a more spherical in-fall. In sim- 
ulation B the minimum moves to higher mass values as the 
merger of the initially three clumps progresses. At the end 
of the simulation the minimum occurs at ~ 1 x 10 5 M©. 



mass of « 1 x 10 5 M© and a radius of « 0.6 pc. 

In figure [10] we show the rotationally supported objects 
face-on and edge-on. We have plotted iso-density surfaces 
using overdensities within a relatively narrow range (approx- 
imately an order of magnitude). These visualisations give a 
clear picture of the shape of the central objects that form 
in all three simulations. Each object is quite clearly a disc. 
The feature in the edge on visualisation of simulation A (top 
panel) which points north-west is the tidal tail seen at the 
bottom of the face-on view. There is another, less promi- 
nent, tidal tail at the top of the disc in simulation A which 
appears at the bottom left in the edge on view. The disc in 
simulation B (middle panel) is less obstructed, the smaller 
clump is also visible in both visualisations. The disc in sim- 
ulation C (bottom panel) is the cleanest example of a disc. 
Like the disc in simulation B it has a somewhat ellipsoidal 
shape. 



4.4 Formation of a self-gravitating massive fat 
disc 

About 0.1-1% of the gas in the inner part of both haloes has 
settled into a rotationally-supported self- gravitating object 
with rotation velocities of 25 — 60kms -1 . We now progress 
to inspect these structures in somewhat more detail. 

In figure [9] we show the final panels of figures [1] and [2] at 
a larger scale. The mass of the central object in simulation 
A in the left hand panel is ~ 2 x 10 4 M©. The radius of the 
central object is ~ 0.3 pc with a compact inner core with 
a radius of ~ 0.1 pc. In the right hand panel we show the 
central object in simulation B at the final output time. The 
mass of the central object (large clump) is « 1 x 10 5 M© 
and the radius is ~ 0.5 pc. The second clump on the right 
hand side has a mass of ~ 5 x 10 3 M© and a radius of ~ 0.1 
pc. The disc in simulation C (shown in figure [3} has a final 



4.5 Properties of the centrifugally supported 
disc(-like object) 

We now have a closer look at the surface mass density profile 
and the level of rotational support of the disc(-like) object 
for each simulation. As can be seen in the left panel of fig- 
ure [TT] in all three cases the surface mass density profile is 
exponential over several scale lengths. The scale lengths are 
0.075, 0.20 and 0.27 pc in simulations A-C, respectively. At 
the outer "edge" of the disc at about 0.3 - 1 pc the surface 
mass density profile reverts to that expected for an isother- 
mal spherical density distribution in which the discs are em- 
bedded. 

In the middle panel of figure [TT] we compare the density 
structure of the discs to that of an exponential disc which is 
given by 
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Figure 10. The density distribution of the rotationally supported fat discs from all three simulations shown 'face-on' (left panels) and 
edge on (right panels). In all cases (approximate) iso-density surfaces are plotted. Note the prominent tidal tail(s) in simulation A and 
the secondary clump in simulation B. The mass of the disc in simulation A is « 2 x 10 4 Mq, while the disc in simulation B has a mass 
of « 1 x 10 5 Mq, with the smaller clump having a mass of ^5 x 1O 3 M . The mass of the disc in simulation C is « 1 x 10 5 Mq. 
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Figure 11. Left-hand Panels: The surface mass density profile of the rotationally supported gas at the centre of each halo shown in 
figure 1101 Middle Panel: The black solid lines are the contours of the azimuthally averaged density in the R — z plane of the simulated 
discs. The red contours and the blue dashed curve are the density and vertical scale height an isothermal exponential disc with the same 
central density, scale length and temperature would have. Right Panel: The four velocities; V ro t, V c , \l\/a\ and \l\/\r\ as a function of 
radius. 



n(r, z) = n exp ^ - ^- J sech 2 y^=^~ J > ( 14 ) 
with scale height 



(47rG/im H n e- 2 -/^d)i/2' ^ ^ 

where c s is the sound speed, \i is the mean molecular weight, 
no is the central density, r is the distance fr om the centre 
of mass and Rg is th e scale-radius of the disc (|Spitzerlll942l ; 
lOh fe Haimemll2002h . We show contours of the azimuthally 
averaged density in the R — z plane of the discs. The discs 
(especially that in simulation A) are somewhat fatter than 
expected if they were supported by pressure only. This is not 
surprising considering the rather unrelaxed dynamical state 
of the discs. The general agreement of the density structure 
with that of an exponential isothermal disc is nevertheless 
remarkable. In simulation B the structure in the top right 
corner is due to the secondary less massive clump. 



In the right panel of figure [TT] we compare the actual ro- 
tation velocities of the gas in the discs (shown by the black 
solid curve) to our two estimates for the rotational veloc- 
ity |/|/|ai| (black dot-dashed curve) and \l\/\r\ (red curve) 
and the circular velocity (blue curve). We compute the ac- 
tual rotational velocities by rotating our coordinate system 
into the coordinate system of the disc using the matrix of 
eigenvectors obtained from the inertia tensor (which we then 
checked visually). The rotation velocity in the plane of the 
disc is then easily calculated using trigonometic transforms. 
The peak of the actual rotation velocities occurs at radii cor- 
responding to, two to three times the scale radius of the ex- 
ponential discs and range from 25 to 60 kms -1 . The actual 
rotation velocities agree with our estimate from the angular 
momentum vector and the largest eigenvalue of the inertia 
tensor, |/|/|ai| to within 10-20%, despite the fact that the 
latter was calculated for the enclosed mass in spheres cen- 
tred on the densest cell. As already mentioned the other 
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estimate for the rotation velocities |/|/|r| is systematically 
lower by a factor 1.5-3. 

The actual rotation velocities exceed the circular veloc- 
ities by up to a factor of 2. This is probably due to a com- 
bination of reasons. The peak rotation velocities of a thin 
exponential disc is about 15 % larger than that of a spher- 
ical mass configuration with the same enclosed mass. More 
importantly the discs in simulation B and C and probably 
also that in simulation A have an ellipsoidal shape which 
should significantly raise the rotation velocities compared 
to the circular velocity of a spherical mass distribution. Fi- 
nally the discs have probably not yet reached centrifugal 
equilibrium. The gas is still settling into rotational support 
and has probably fallen to somewhat smaller radius initially 
than expected for centrifugal equilibrium and will take some 
time to reach centrifugal equilibrium. 

4.6 Numerical Limitations 

Probing the collapse of gas at the centre of dark matter 
haloes in a cosmological context is an extremely challenging 
exercise when we wish to follow the collapse to very high 
densities. In this work we have reached radii as small as 
0.01% of the virial radius while at the same time following 
the dynamical evolution of a substantial fraction of the gas 
in the halo. This presents a considerable challenge even for 
an AMR code like Enzo. The main problem is that the high 
refinement levels necessary to achieve such a large dynamic 
range mean that the code will normally grind to a halt fol- 
lowing the dynamical evolution of whatever is the first high 
density region to form. As discussed, in order to ameliorate 
this problem we have changed the refinement level during 
the simulation in order to allow our three haloes of choice 
to build up to their full mass. To what level may this have 
affected our results? 

In order to address this we have run simulations with 
twice the initial refinement level and sixty four times the 
initial refinement level and found no systematic difference 
in the initial value of the angular momentum of the halo. 
While the detailed dynamical evolution is obviously differ- 
ent especially in the later stages of the collapse when the 
initial distribution of matter in the halo is different the qual- 
itative behavior of the dynamical evolution did not change. 
We would also like to point out that our results are in most 
aspects similar to those of W07 who did not change the re- 
finement level during the simulation. This meant, however, 
that they were not able to follow the gas at the centre of the 
halo settling into rotational support. We feel that this is the 
most interesting aspect of our work and well worth explor- 
ing the limits of the code with some manual intervention. 
Obviously further work is needed to address these questions 
in more detail. 



5 CONCLUSIONS 

We have investigated the dynamical evolution of the gas in 
haloes with virial temperatures of between ~ 13000 K and 
~ 30000 K assuming that cooling by atomic hydrogen and 
helium are the dominant cooling processes. The dynamical 
evolution of the gas in the haloes is complex and highly tur- 
bulent with turbulent velocities approaching the virial ve- 



locity of the haloes. The gas in the inner part of our haloes 
collapses isothermally with a temperature of 6000-7000K on 
a somewhat slower time scale than the free-fall time and 
settles into a close to isothermal (p oc r~ 2 ) density profile. 
We find no signs of efficient fragmentation confirming sug- 
gestions that the isothermal collapse of gas in a dark matter 
halo at temperatures close to the virial temperature of the 
halo leads to only modest fragmentation. The inner 0.1-1 
% of the gas in the virialised haloes loses as much as 95% 
of its initial angular momentum during the collapse. The 
gas thereby collapses by a factor of 300-1000 in radius and 
eventually settles into a very compact rotationally supported 
self- gravitating disc with peak rotation velocities of 25 - 60 
kms" 1 and "radii" of 0.3 pc - 0.6 pc (0.05% - 0.1% of the 
virial radius of the host halo). 

The discs have an exponential surface mass density pro- 
file with scale length in the range 0.075 — 0.27 pc which ex- 
tends over several scale lengths. The vertical structure of 
the disc is somewhat more extended than expected for a 
purely pressure supported isothermal axisymmetric expo- 
nential disc. 

Massive compact self- gravitating discs such as those 
found in our simulations have been suggested to evolve 
into massive seed black holes which later in the hierarchical 
build-up of galaxies will grow into the supermassive black 
holes found at the centre of the bulges of present-day galax- 
ies. Unfortunately we were not yet able to follow the further 
dynamical evolution of the discs or the gas further out in the 
haloes in our simulations. However, independent of whether 
the gas in these discs will continue to efficiently lose angu- 
lar momentum and contract further or will fragment and 
form an ultra- compact star cluster we most probably have 
identified an important intermediary stage en route to the 
formation of a massive seed black hole. 
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